/*!  \file  MRRR.cpp
 *	  \brief  
 * 
 *   \details   details
 *   \author    Yordan Madzhunkov German Research School for Simulation Sciences GmbH
 *   \version   1.0
 *   \date      06/29/2013 11:26:33 PM
 *   \copyright GNU Public License.
 */
#include <stdlib.h> // setenv
#include <stdio.h> // sprintf
#include "MRRR.h"
#include "Matrix.h"
extern "C" {
/* Parallel computation of all or a subset of eigenvalues and 
 * optionally eigenvectors of a symmetric tridiagonal matrix based on 
 * the algorithm of Multiple Relatively Robust Representations (MRRR). 
 * The routine targets multi-core architectures. 
 * The implementation is based on LAPACK's routine 'dstemr'.
 *
 * Function prototype: */

int mrrr(char *jobz, char *range, int *n, double * d,
	 double * e, double *vl, double *vu, int *il, int *iu,
	 int *tryrac, int *m, double *W, double *Z, int *ldz,
	 int *Zsupp);

/* Arguments:
 * ----------
 *
 * INPUTS: 
 * -------
 * jobz              "N" or "n" - compute only eigenvalues
 *                   "V" or "v" - compute also eigenvectors
 *                   "C" or "c" - count the maximal number of 
 *                                locally computed eigenvectors
 * range             "A" or "a" - all
 *                   "V" or "v" - by interval: (VL,VU]
 *                   "I" or "i" - by index:     IL-IU
 * n                 Matrix size
 * ldz               Leading dimension of eigenvector matrix Z; 
 *                   often equal to matrix size n
 *
 * INPUT + OUTPUT: 
 * ---------------
 * D (double[n])     Diagonal elements of tridiagonal T.
 *                   (On output the array will be overwritten).
 * E (double[n])     Off-diagonal elements of tridiagonal T.
 *                   First n-1 elements contain off-diagonals,
 *                   the last element can have an abitrary value. 
 *                   (On output the array will be overwritten.)
 * vl                If range="V", lower bound of interval
 *                   (vl,vu], on output refined.
 *                   If range="A" or "I" not referenced as input.
 *                   On output the interval (vl,vu] contains ALL
 *                   the computed eigenvalues.
 * vu                If range="V", upper bound of interval
 *                   (vl,vu], on output refined.
 *                   If range="A" or "I" not referenced as input.
 *                   On output the interval (vl,vu] contains ALL
 *                   the computed eigenvalues.
 * il                If range="I", lower index (1-based indexing) of 
 *                   the subset 'il' to 'iu'.
 *                   If range="A" or "V" not referenced as input.
 *                   On output the eigenvalues with index il to iu are 
 *                   computed.
 * iu                If range="I", upper index (1-based indexing) of 
 *                   the subset 'il' to 'iu'.
 *                   If range="A" or "V" not referenced as input.
 *                   On output the eigenvalues with index il to iu are 
 *                   computed.
 * tryrac            0 - do not try to achieve high relative accuracy.
 *                   1 - relative accuracy will be attempted; 
 *                       on output it is set to zero if high relative 
 *                       accuracy is not achieved.
 *
 * OUTPUT: 
 * -------
 * m                 Number of eigenvalues and eigenvectors computed. 
 *                   If jobz="C", 'm' will be set to the number of 
 *                   eigenvalues/eigenvectors that will be computed. 
 * W (double[n])     Eigenvalues
 *                   The first 'm' entries contain the eigenvalues.
 * Z                 Eigenvectors.
 * (double[n*m])     Enough space must be provided to store the
 *                   vectors. 'm' should be bigger or equal 
 *                   to 'n' for range="A" and 'iu-il+1' for range="I".
 *                   For range="V" the minimum 'm' can be queried using
 *                   jobz="C".
 * Zsupp             Support of eigenvectors, which is given by
 * (double[2*n])     Z[2*i] to Z[2*i+1] for the i-th eigenvector
 *                   stored locally (1-based indexing).
 *
 * RETURN VALUE: 
 * -------------
 *                 0 - Success  
 *                 1 - Wrong input parameter
 *                 2 - Misc errors  
 *
 *
 */
}

MRRR::MRRR() : TriEigen("Parallel mrrr") {
	setCompileTime(__DATE__, __TIME__); 
	canComputeRange = true;
	canBeParallel = true;
}
bool MRRR::setNumberOfThreads(const int nThreads_){ 
	char buf[16]; sprintf(buf, "%d",  nThreads_);
	const int overwrite = 1;
	if (setenv("PMR_NUM_THREADS", buf, overwrite) == 0) {
		nThreads = nThreads_;
		cout << "PMR_NUM_THREADS = " << getenv("PMR_NUM_THREADS") << endl;
		return true;
	} else 
		return false;
};
void MRRR::compute(const Matrix & T, Matrix & Sigma, Matrix & Z ){
	reset();
	assert(T.isTriDiagonal() && Sigma.isDiagonal() && Z.isOrthogonal());
	char jobz = 'V', range = 'A';
	int n = T.cols(), info = -1, ldz = Z.lda();
	double vl = 0, vu = 1; int il = 0, iu = n-1;
	int m = 0; // total number of eidgen values found
	double * w = Sigma.getModData();
	double * z = Z.getModData(); 
	if (Z.isRowMajor()) Z.makeColMajor();
	double * d = NULL, *e = NULL; allocate(d, n); allocate(e, n);
	memcpy(d, T.getData()       ,    n * sizeof(*d));
	memcpy(e, T.getOffDiagonal(), (n-1)* sizeof(*d));
	e[n-1] = 0;	int tryrac = 1;
	int * isupp = NULL; allocate(isupp, 2 * n); 
	info = mrrr(&jobz, &range, &n, d, e, &vl, &vu, &il, &iu, &tryrac, &m, 
			w, z, &ldz, isupp);
	success &= info == 0;
	delocate(isupp); delocate(e); delocate(d);
};
void MRRR::computeRange( const Matrix & T, Matrix & Sigma, Matrix & Z, 
	const int indexLow, const int indexUp) {
	reset();
	assert(T.isTriDiagonal() && Sigma.isDiagonal() && Z.isOrthogonal());
	char jobz = 'V', range = 'I';
	int n = T.cols(), info = -1, ldz = Z.lda();
	double vl = 0, vu = 1; int il = indexLow, iu = indexUp;
	int m = 0; // total number of eidgen values found
	double * w = NULL; allocate(w, n);//Sigma.getModData();
	double * z = Z.getModData(); 
	if (Z.isRowMajor()) Z.makeColMajor();
	double * d = NULL, *e = NULL; allocate(d, n); allocate(e, n);
	memcpy(d, T.getData()       ,    n * sizeof(*d));
	memcpy(e, T.getOffDiagonal(), (n-1)* sizeof(*d));
	e[n-1] = 0;	int tryrac = 1;
	int * isupp = NULL; allocate(isupp, 2 * n); 
	info = mrrr(&jobz, &range, &n, d, e, &vl, &vu, &il, &iu, &tryrac, &m, 
			w, z, &ldz, isupp);
	memcpy(Sigma.getModData(), w, (iu - il + 1) * sizeof(*w));
	success &= info == 0;
	delocate(w);
	delocate(isupp); 
	delocate(e); 
	delocate(d);
}

BiNEMRRR::BiNEMRRR() : BiSVD( "NE Parallel MRRR" ) {
	aUpper = true; 
	aLower = false; 
	aVColMajor = true; 
	aUColMajor = true; 
	VColVec = false;
	setTriEigen(*(new MRRR()));
	setCompileTime(__DATE__, __TIME__); 
	if (eigen) setCompileTime(eigen->getCompileTime());
};
BiNE2MRRR::BiNE2MRRR() : BiSVD( "NE2 Parallel MRRR" ) {
	aUpper = true; 
	aLower = false; 
	aVColMajor = true; 
	aUColMajor = true; 
	VColVec = false;
	setTriEigen(*(new MRRR()));
	setCompileTime(__DATE__, __TIME__);
	if (eigen) setCompileTime(eigen->getCompileTime());
};
BiGKMRRR::BiGKMRRR() : BiSVD( "GK Parallel MRRR" )  {
	aUpper = true; 
	aLower = false; 
	aVColMajor = true; 
	aUColMajor = true; 
	VColVec = false;
	setTriEigen(*(new MRRR()));
	setCompileTime(__DATE__, __TIME__);
	if (eigen) setCompileTime(eigen->getCompileTime());
}
REGISTER(MRRR);
REGISTER(BiNEMRRR );
REGISTER(BiNE2MRRR );
REGISTER(BiGKMRRR );
